Data Analysis and Visualization
Paper 2: Community Data for 4 countries
WORK IN PROGRESS

Section 1. Data cleaning and Demographics

Display code

options(digits = 2)
rm(list=ls())

## Install "pacman" package if not installed
# (remove the # symbol from the line below):
# install.packages("pacman")

## Load R packages:
pacman::p_load(data.table, tidyverse, haven, labelled, vtable, 
               psych, scales, weights, clipr, forcats,
               stargazer, ggthemes, ggcharts, geomtextpath,
               corrplot, tm, readxl, patchwork, modelsummary,
               gt, lme4, car, lmerTest, 
               ggeffects, magrittr, broom, broom.mixed,
               backports, effects, interactions, plyr, sjPlot)

## Import datasets:
## These are newer datasets with new variables created in the individual analyses:
dsgmb <- fread("~/Dropbox/2023_Gagan/Paper2/data/dsgmb.csv")
dsgmb$Country <- "Gambia"
dsgmb$id <- 1:nrow(dsgmb)
dsgmb$ID <- paste0("GMB",dsgmb$id)

dspak <- fread("~/Dropbox/2023_Gagan/Paper2/data/dspak.csv")
dspak$Country <- "Pakistan"
dspak$id <- 1:nrow(dspak)
dspak$ID <- paste0("PAK",dspak$id)

dstza <- fread("~/Dropbox/2023_Gagan/Paper2/data/dstza.csv")
dstza$Country <- "Tanzania"
dstza$id <- 1:nrow(dstza)
dstza$ID <- paste0("TZA",dstza$id)

dsuga <- fread("~/Dropbox/2023_Gagan/Paper2/data/dsuga.csv")
dsuga$Country <- "Uganda"
dsuga$id <- 1:nrow(dsuga)
dsuga$ID <- paste0("UGA",dsuga$id)

## Correct asterisk pattern for stargazer tables:
starpattern <- "<em>&#42;p&lt;0.05;&nbsp;&#42;&#42;p&lt;0.01;&nbsp;&#42;&#42;&#42;p&lt;0.001</em>"

## List of variables to retain from all datasets:
list1 <- c("ID", "Country", "age", "gender", 
           "ses", "jobnature", "religion", "married", "education", "Ethnicity",
           "IGF01", "IGF02", "IGF03", "IGI01", "IGI02", "IGI03", 
           "OGF01", "OGF02", "OGF03", "OGI01", "OGI02", "OGI03",
           "ENDBCL01", "ENDBCL02", "ENDBCL03", "ENDBBL01", "ENDBBL02", "ENDBBL03",
           "EXPBCL01", "EXPBCL02", "EXPBCL03", "EXPBBL01", "EXPBBL02", "EXPBBL03", 
           "empathic_concern_01", "empathic_concern_02", "empathic_concern_03",
           "perspective_taking_01", "perspective_taking_02", 
           "perspective_taking_03", "perspective_taking_04", "history_discrimination",
           "og_hostility", "og_cooperation", "fight_outgroup", "imagistic",
           "event_positive_affect", "event_negative_affect", "event_episodic_recall",
           "event_shared_perception", "event_event_reflection", 
           "event_transformative_indiv", "event_transformative_group",
           "freq_positive_contact", "freq_negative_contact", 
           "sprf", "exp_religious_freedom")

## Subset datasets to only the columns in the above list:
dsgmb1 <- dsgmb[, ..list1]
dspak1 <- dspak[, ..list1]
dstza1 <- dstza[, ..list1]
dsuga1 <- dsuga[, ..list1]

## Merged dataset with needed columns only
ds <- rbind(dsgmb1, dspak1, dstza1, dsuga1)

## Rename the "Event_" columns with title case:
ds01 <- ds[, 1:44]
ds02 <- ds[, !(2:44)]
colnames(ds02) <- stringr::str_to_title(colnames(ds02))
ds02$ID <- ds02$Id
ds02 <- ds02[, !1]
ds <- merge(ds01, ds02, by = "ID")
rm(ds01, ds02)

Data cleaning

Criteria:

• Survey duration:

not available for this dataset as it was manual interviews

• Bad open ended responses:

not applicable since all are manual responses

• SD on unrelated scales

• INCOME VARIABLE

some cases with outlier responses in “income” variable need to be removed. THIS PART NEEDS TO BE REDONE

• Missing data in scales

Missing data:

Display code

mds <- map(ds, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       2       2       3       2      34 

Display code

mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data")+
  theme_bw()

lp01

Display code

lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data")+
  theme_bw()

lp02

Display code

## Individual country datasets:
dsgmb <- ds[ds$Country == "Gambia", ]
dspak <- ds[ds$Country == "Pakistan", ]
dstza <- ds[ds$Country == "Tanzania", ]
dsuga <- ds[ds$Country == "Uganda", ]

Missing data: Gambia

Display code

mds <- map(dsgmb, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.4     0.8     1.3     1.7     9.3 

Display code

mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Gambia")+
  theme_bw()

lp01

Display code

lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Gambia")+
  theme_bw()

lp02

Missing data: Pakistan

Display code

mds <- map(dspak, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       2       0     100 

Display code

mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp01

Display code

lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp02

Pak after removing BCL01:

Display code

dspak22 <- subset(dspak, select = -c(EXPBCL01))

mds <- map(dspak22, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.00    0.00    0.07    0.00    0.79 

Display code

mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:54,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp01

Display code

lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp02

Missing data: Tanzania

Display code

mds <- map(dstza, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       5       6       7       8      44 

Display code

mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Tanzania")+
  theme_bw()

lp01

Display code

lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Tanzania")+
  theme_bw()

lp02

Missing data: Uganda

Display code

mds <- map(dsuga, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.7     0.9     1.1     1.3     5.1 

Display code

mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Uganda")+
  theme_bw()

lp01

Display code

lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Uganda")+
  theme_bw()

lp02

Demographic Variables:

Display code

ds$gender <- ifelse(ds$gender == "Male", "Male",
             ifelse(ds$gender == "Female", "Female", NA))

## ses:
ds$ses <- ifelse(ds$ses == "", NA, ds$ses)

## jobnature:
ds$jobnature <- ifelse(ds$jobnature == "", NA, ds$jobnature)

#sentence case for jobnature:
ds$jobnature <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$jobnature, perl = TRUE)

ds$jobnature <- ifelse(ds$jobnature == "Non-government/self-employed", 
                       "Non-government", ds$jobnature)

## Religion:
ds$religion <- ifelse(ds$religion == "", NA, ds$religion)

ds$religion <- ifelse(ds$religion == "Christian (Catholic)", "Christian: Catholic", 
               ifelse(ds$religion == "Christian (Protestant)", "Christian: Protestant",
               ifelse(ds$religion == "Muslim (Shia)", "Muslim: Shia",
               ifelse(ds$religion == "Muslim (Sunni)", "Muslim: Sunni", ds$religion))))

## Marital status

ds$married <- ifelse(ds$married == "Not married", "Unmarried", ds$married)
ds$married <- ifelse(ds$married == "", NA, ds$married)

## Individual country datasets:
dsgmb <- ds[ds$Country == "Gambia", ]
dspak <- ds[ds$Country == "Pakistan", ]
dstza <- ds[ds$Country == "Tanzania", ]
dsuga <- ds[ds$Country == "Uganda", ]

Demographic variables in the analysis:

• Age
• Gender
• Socio-economic status
• Nature of employment
• Religious affiliation
• Marital status
• Ethnicity

Variable: Sample size by Country

Display code

tbl01 <- table(ds$Country)

## Table of user language by country:
tbl01

  Gambia Pakistan Tanzania   Uganda 
     237      505      352      448 

Display code

## Sample size by country:

lp01 <- ds %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = Country,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Sample size by country")+
  theme_bw()

lp01

Variable: Age

Display code

summary(ds$age)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
     18      26      33      37      45      92      29 

Display code

ds %>% drop_na(age)%>%
ggplot(aes(x = age))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 37.00", 
                 xintercept = 37.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Age", 
       y = "Frequency", 
       title = "Age distribution (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(age)%>%
ggplot(aes(x = age))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  labs(x = "Age", 
       y = "Frequency", 
       title = "Age distribution by country")+
  facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Gender

Display code

ds$gender <- ifelse(ds$gender == "Male", "Male",
             ifelse(ds$gender == "Female", "Female", NA))

lp02 <- ds %>% 
  drop_na(gender, age, Country) %>%
lollipop_chart(x = gender,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Gender distribution (full sample)")+
#  facet_wrap(~Country, nrow = 2)+
  theme_bw()

lp02

Display code

bp01 <- ds %>% drop_na(gender, age) %>% 
  ggplot(aes(y = age, 
             x = gender))+
geom_boxplot(fill = "grey")+
  labs(y = "Age",
       x = "",
       title = "Age and gender distribution by country")+
  facet_wrap(~Country, nrow = 2)+
  coord_flip()+
  theme_bw()

bp01

Variable: Socio-economic status

Display code

ds$ses <- ifelse(ds$ses == "", NA, ds$ses)
table(ds$ses)

Lower middle       Middle        Upper Upper middle 
         479          831           27          172 

Display code

ds %>% drop_na(ses) %>%
lollipop_chart(x = ses,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Socioeconomic status (full sample)")+
  theme_bw()

Variable: Nature of employment

Display code

table(ds$jobnature)

    Government Non-government          Other        Retired  Self-employed 
           329            295              2              1            773 
       Student 
             6 

Display code

ds %>% drop_na(jobnature) %>%
lollipop_chart(x = jobnature,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Nature of employment (full sample)")+
  theme_bw()

Variable: Religious affiliation

Display code

table(ds$religion)

    Christian: Catholic        Christian: Other  Christian: Pentecostal 
                    460                      43                      33 
  Christian: Protestant                   Hindu           Muslim: Other 
                    209                      37                       4 
           Muslim: Shia           Muslim: Sunni                   Other 
                    188                     433                      87 
                   Sikh Spiritual not religious 
                      9                      14 

Display code

lp05 <- ds %>% drop_na(religion) %>%
lollipop_chart(x = religion,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Religious distribution (full sample)")+
  theme_bw()

lp05

Variable: Marital status

Display code

table(ds$married)

  Married     Other Unmarried 
      957        94       483 

Variable: Education

Display code

table(ds$education)

                                         
                                      44 
                                       4 
                                       8 
                                   4m  4 
                                       1 
                                       7 
                                       7 
                         Advance Diploma 
                                       1 
     Advance Diploma in Business Studies 
                                       1 
           Advance Diploma in Management 
                                       3 
              Advance Diploma in Nursing 
                                       1 
   Advance Diploma in Tourism Management 
                                       1 
                           Advance level 
                                       1 
                                   B.S.C 
                                       1 
                                Bachelor 
                                     104 
                         Bachelor Degree 
                                       2 
     Bachelor in Business Administration 
                                       1 
                       Bachelor's Degree 
                                       3 
  Bachelor's Degree in Political Science 
                                       1 
                     Bachelors (ongoing) 
                                       1 
                        Bachelors Degree 
                                       1 
                  Bachelors in Education 
                                       1 
                       Bachelors of Edu. 
                                       1 
                                     Bsc 
                                       1 
                                     BSc 
                                       3 
                                     BSC 
                                       5 
                  BSc (Arabic Education) 
                                       1 
                   BSc (Hons) Management 
                                       1 
                       BSc in Accounting 
                                       1 
                         BSc in Economic 
                                       1 
                        BSc in Economica 
                                       1 
                        BSc in Education 
                                       1 
                   BSc in Human Resource 
                                       1 
                       BSc in Management 
                                       1 
                BSc in Political Science 
                                       1 
                      BsSc of Accounting 
                                       1 
                             Certificate 
                                       3 
   Certificate in International Relation 
                                       1 
       Certificate in Management Studies 
                                       1 
                       Certificate Level 
                                       1 
    Certificate on Animal Health Product 
                                       1 
                            Certtificate 
                                       1 
                        Chartered Banker 
                                       1 
                                   Cheti 
                                       6 
                                    Chuo 
                                       6 
                             Chuo Ckikuu 
                                       1 
                              Chuo kikuu 
                                      11 
                              Chuo Kikuu 
                                       3 
                             Claiasa nne 
                                       1 
                                 College 
                                       9 
                           College (HTC) 
                                       1 
                           College (HWD) 
                                       1 
                     College Certificate 
                                       1 
                         currently in s4 
                                       1 
                                Darasa 7 
                                       8 
                                Darasa 8 
                                       1 
                             Darasa la 7 
                                       1 
                          Darasa la Saba 
                                       1 
                              Darsa la 7 
                                       2 
                              Darsa La 7 
                                       3 
                           Darsa la saba 
                                       1 
                                  Degree 
                                      15 
               Degree in Islamic Studies 
                                       1 
                          Degree-Shahada 
                                       1 
                           Didato channe 
                                       1 
                           Digri/Shahada 
                                       1 
                                 Diploma 
                                      22 
       Diploma 2 Business Administration 
                                       1 
           Diploma in Community Policing 
                                       1 
                          Diploma in ICT 
                                       2 
                   Diploma in Management 
                                       3 
                     Diploma in Theology 
                                       1 
           Diploma in Travel and Tourism 
                                       1 
                                Diplomat 
                                       1 
                            Elimu Msingi 
                                       1 
           Form 3 (Old education System) 
                                       2 
           Form 3 (Old Education System) 
                                       1 
                                  Form 4 
                                      15 
                                  Form 5 
                                       1 
           Form 5 (Old education System) 
                                       2 
           Form 5 (Old Education System) 
                                       1 
                                  Form 6 
                                       3 
                               Form Five 
                                       2 
                             Form Five 5 
                                       2 
                               Form Four 
                                      73 
                             Form Four 4 
                                       3 
                             Form Four 5 
                                       1 
                                 Form IV 
                                       1 
                              Form One 1 
                                       1 
                                Form Six 
                                      14 
                              Form Three 
                                       1 
                                Form Two 
                                       6 
                              Form Two 2 
                                       1 
                    Gambia College (PTC) 
                                       1 
                       GCE Level(WASSCE) 
                                       1 
                  Grade 8 ( Upper Basic) 
                                       1 
                        Graduate Diploma 
                                       1 
                             High School 
                                       4 
                    High School (WASSCE) 
                                       1 
                 High School Certificate 
                                       2 
                       High School Level 
                                       1 
           High School-Islamic Education 
                                       1 
                          Higher Diploma 
                                       1 
       Higher Diploma in Education (HDC) 
                                       1 
        Higher Teacher Certificate (HTC) 
                                       7 
      Higher Teacher's Certificate (HTC) 
                                       5 
                            Hotel School 
                                       1 
                                     HTC 
                                       2 
                             HTC Primary 
                                       1 
                                 HTC/HNC 
                                       1 
                              Illiterate 
                                      10 
            Informal Education (Madrasa) 
                                       1 
                            Intermediate 
                                      85 
                       Isachetovs Degree 
                                       1 
                          Islamic School 
                                       1 
              Islamic Traditional School 
                                       1 
                                      IV 
                                       1 
                       Junior Sec School 
                                       1 
                        Junior Secondary 
                                       1 
                                Kidato 4 
                                       3 
                            Kidato cha 1 
                                       1 
                            Kidato cha 4 
                                       1 
                            Kidato cha 6 
                                       2 
                          Kidato cha nne 
                                       2 
                         Kidato Cha seta 
                                       1 
                           Kidato Channe 
                                       1 
                           Kidato chaune 
                                       1 
                                    La 4 
                                       3 
                                    La 7 
                                       7 
                                    LA 7 
                                       8 
                                  La nne 
                                       4 
                              La nne (4) 
                                       1 
                          La nne mkoloni 
                                       1 
                                 La saba 
                                      17 
                               La saba 7 
                                       8 
                           Lanne Mkoloni 
                                       2 
                                  Lasaba 
                                       1 
                                     LLB 
                                       1 
                                  Master 
                                     161 
                          Masters Degree 
                                       2 
                    Masters in Economics 
                                       1 
                    Masters in Sociology 
                                       1 
                             Matriculate 
                                      34 
Memorized the Quran (Informal Education) 
                                       1 
                                  Middle 
                                      53 
                                   MPhil 
                                      56 
                                     MSC 
                                       1 
                                  Msingi 
                                       3 
                                    None 
                                      50 
                                 O Level 
                                       3 
                         Ordinary Levels 
                                       1 
                                     PhD 
                                       3 
                           Post Graduate 
                                       2 
                   Post Graduate Diploma 
                                       1 
                                 Primary 
                                     213 
                          Primary School 
                                       3 
       Primary Teacher Certificate (PTC) 
                                       2 
     Primary Teacher's Certificate (PTC) 
                                       1 
                                     PTC 
                                       1 
            Quranic High School Graduate 
                                       1 
                          Quranic School 
                                       1 
                              Seconary 4 
                                       1 
                               Secondary 
                                     143 
                         Secondary  Four 
                                       1 
                         Secondary level 
                                       1 
                        Secondary School 
                                       6 
                                Sekonorw 
                                       1 
                Senior School Graduation 
                                       1 
                        Senior Secondary 
                                       2 
                 Senior Secondary School 
                                      16 
    Senior Secondary(WASSCE Certificate) 
                                       1 
                                 Shadada 
                                       3 
                                 Shahada 
                                       5 
                             Shahada yai 
                                       1 
                             Shule Msiry 
                                       1 
                         Shule ya msingi 
                                       1 
                                Sokondal 
                                       1 
                              StaShahada 
                                       2 
                                 Std VII 
                                       3 
                                 STD VII 
                                       8 
                                Tertiary 
                                      33 
                      Tertiary Education 
                                      10 
                          Tertiary level 
                                       1 
             Traditional Quranic Schools 
                                       1 
                            Udergraduate 
                                       1 
  Undergraduate Degree in Communications 
                                       1 
                              University 
                                      15 
                       University Degree 
                                       7 
                University of the Gambia 
                                       1 
                                 Uzamili 
                                       1 
                                      VI 
                                       1 
                                 VII Std 
                                       1 
                              Vocational 
                                       8 
                                  WASSCE 
                                      16 
          WASSCE examination certificate 
                                       1 

Variable: Ethnicity

Display code

## Gambia:
eth <- as.data.frame(table(dsgmb$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 2, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dsgmb2 <- dsgmb[, c("Ethnicity")]

dsgmb2$ethnicity <- ifelse(dsgmb2$Ethnicity %in% l1, dsgmb2$Ethnicity, "Other")
dsgmb2$ethnicity <- ifelse(dsgmb2$ethnicity == "Serere", "Serer",
                 ifelse(dsgmb2$ethnicity == "", NA, dsgmb2$ethnicity))

table(dsgmb2$ethnicity)

 Balanta     Fula     Jola Karonika Mandinka  Manjago    Other Sarahula 
       3       36       46        3       88        5       16        2 
   Serer   Wollof 
       5       25 

Display code

ethgmb <- dsgmb2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Gambia")+
  theme_bw()

ethgmb

Display code

## Pakistan:
eth <- as.data.frame(table(dspak$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 7, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dspak2 <- dspak[, c("Ethnicity")]

dspak2$ethnicity <- ifelse(dspak2$Ethnicity %in% l1, dspak2$Ethnicity, "Other")
dspak2$ethnicity <- ifelse(dspak2$ethnicity == "Serere", "Serer",
                 ifelse(dspak2$ethnicity == "", NA, dspak2$ethnicity))

table(dspak2$ethnicity)

      Bhatti    Christain       Hazara         Jopu        Other      Punjabi 
          22           28            7           20          124          118 
      Rajput      Sadozai         Shia         Sikh       Sindhi        Sunni 
           9           10           20           10           14           29 
        Syed Urduspeaking 
          34           60 

Display code

ethpak <- dspak2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Pakistan")+
  theme_bw()

ethpak

Display code

## Tanzania:
eth <- as.data.frame(table(dstza$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 7, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dstza2 <- dstza[, c("Ethnicity")]

dstza2$ethnicity <- ifelse(dstza2$Ethnicity %in% l1, dstza2$Ethnicity, "Other")
dstza2$ethnicity <- ifelse(dstza2$ethnicity == "Serere", "Serer",
                 ifelse(dstza2$ethnicity == "", NA, dstza2$ethnicity))

table(dstza2$ethnicity)

 Maasai  Mchaga   Mgogo Mlugulu Mluguru  Mngoni   Mpare Msambaa Msukuma    Muha 
    122       9       8      11      15       8       9      11      10       7 
 Mzigua   Other 
      7     120 

Display code

ethtza <- dstza2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Tanzania")+
  theme_bw()

ethtza

Display code

## Uganda:
eth <- as.data.frame(table(dsuga$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 7, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dsuga2 <- dsuga[, c("Ethnicity")]

dsuga2$ethnicity <- ifelse(dsuga2$Ethnicity %in% l1, dsuga2$Ethnicity, "Other")
dsuga2$ethnicity <- ifelse(dsuga2$ethnicity == "Serere", "Serer",
                 ifelse(dsuga2$ethnicity == "", NA, dsuga2$ethnicity))

table(dsuga2$ethnicity)

     Adhola    Amudaama      Atesot       Bantu      Etesot      Itesot 
         18          11          61          11          56          31 
  Japadhola japhadollah Japhadollah     Munyole      Musoga     Nilotic 
         97          22          23           7           8          27 
      Other 
         76 

Display code

ethuga <- dsuga2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Uganda")+
  theme_bw()

ethuga

Section 2. Measures of interest

Variables of interest in the analysis:

• Endorsement of Barrier Bound Leadership (BBL)

• Endorsement of Barrier Crossing Leadership (BCL)

• Ingroup fusion

• Ingroup identification

• Outgroup bonds (sum of outgroup fusion+identification)

• Empathic concern

• Perspective taking

• Perceived history of discrimination

• Positive contact with outgroup

• Negative contact with outgroup

• Imagistic items:

  • Event: Positive affect
  • Event: Negative affect
  • Event: Episodic recall
  • Event: Shared perception
  • Event: Reflection
  • Event: Transformative for individual
  • Event: Transformative for group
  • Imagistic event (sum of all above)

• Perception of religion Freedom

• Experience of religious freedom

Variable: Endorsement of Barrier Bound Leadership (BBL)

Display code

ds$bbl <- (ds$ENDBBL01+ ds$ENDBBL02+ ds$ENDBBL03)/3
ds$bcl <- (ds$ENDBCL01+ ds$ENDBCL02+ ds$ENDBCL03)/3

summary(ds$bbl)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       3       5       4       6       7      42 

Display code

ds %>% drop_na(bbl)%>%
ggplot(aes(x = bbl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 4.00", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Endorsement of BBL score", 
       y = "Frequency", 
       title = "Endorsement of BBL (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(bbl, Country)%>%
ggplot(aes(x = bbl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Endorsement of BBL score", 
       y = "Frequency", 
       title = "Endorsement of BBL by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(bbl)%>%
ggplot(aes(x = bbl,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Endorsement of BBL score", 
       y = "Frequency", 
       title = "Endorsement of BBL by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Endorsement of Barrier Crossing Leadership (BCL)

Display code

summary(ds$bcl)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      38 

Display code

ds %>% drop_na(bcl)%>%
ggplot(aes(x = bcl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Endorsement of BCL score", 
       y = "Frequency", 
       title = "Endorsement of BCL (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(bcl, Country)%>%
ggplot(aes(x = bcl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Endorsement of BCL score", 
       y = "Frequency", 
       title = "Endorsement of BCL by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(bcl)%>%
ggplot(aes(x = bcl,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Endorsement of BCL score", 
       y = "Frequency", 
       title = "Endorsement of BCL by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Ingroup fusion

Display code

ds$igfusion <- (ds$IGF01+ds$IGF02+ds$IGF03)/3

summary(ds$igfusion)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       6       6       6       7       7      48 

Display code

ds %>% drop_na(igfusion)%>%
ggplot(aes(x = igfusion))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup fusion score", 
       y = "Frequency", 
       title = "Ingroup fusion (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(igfusion, Country)%>%
ggplot(aes(x = igfusion))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Ingroup fusion score", 
       y = "Frequency", 
       title = "Ingroup fusion by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(igfusion)%>%
ggplot(aes(x = igfusion,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup fusion score", 
       y = "Frequency", 
       title = "Ingroup fusion by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Ingroup identification

Display code

ds$IGI01 <- ifelse(ds$IGI01 %in% 1:7, ds$IGI01, NA)
ds$igidentification <- (ds$IGI01+ds$IGI02+ds$IGI03)/3

summary(ds$igidentification)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       6       6       6       7       7      48 

Display code

ds %>% drop_na(igidentification)%>%
ggplot(aes(x = igidentification))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup identification score", 
       y = "Frequency", 
       title = "Ingroup identification (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(igidentification, Country)%>%
ggplot(aes(x = igidentification))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Ingroup identification score", 
       y = "Frequency", 
       title = "Ingroup identification by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(igidentification)%>%
ggplot(aes(x = igidentification,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup identification score", 
       y = "Frequency", 
       title = "Ingroup identification by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Outgroup bonds

Display code

ds$ogbonds <- (ds$OGI01+ds$OGI02+ds$OGI03+
               ds$OGF01+ds$OGF02+ds$OGF03)/6

summary(ds$ogbonds)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       2       3       3       4       7     102 

Display code

ds %>% drop_na(ogbonds)%>%
ggplot(aes(x = ogbonds))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 35)+
  geom_textvline(label = "Mean = 3.00", 
                 xintercept = 3.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Outgroup bonds score", 
       y = "Frequency", 
       title = "Outgroup bonds (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(ogbonds, Country)%>%
ggplot(aes(x = ogbonds))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Outgroup bonds score", 
       y = "Frequency", 
       title = "Outgroup bonds by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(ogbonds)%>%
ggplot(aes(x = ogbonds,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 3.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Outgroup bonds score", 
       y = "Frequency", 
       title = "Outgroup bonds by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Empathic concern

Display code

ds$empathic_concern <- (ds$empathic_concern_01+ds$empathic_concern_02+
                        ds$empathic_concern_03)/3

summary(ds$empathic_concern)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      33 

Display code

ds %>% drop_na(empathic_concern)%>%
ggplot(aes(x = empathic_concern))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 35)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(empathic_concern, Country)%>%
ggplot(aes(x = empathic_concern))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(empathic_concern)%>%
ggplot(aes(x = empathic_concern,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Perspective taking

Display code

ds$perspective_taking <- (ds$perspective_taking_01+ds$perspective_taking_02+
                          ds$perspective_taking_03+ds$perspective_taking_04)/4

summary(ds$perspective_taking)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      54 

Display code

ds %>% drop_na(perspective_taking)%>%
ggplot(aes(x = perspective_taking))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 35)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(perspective_taking, Country)%>%
ggplot(aes(x = perspective_taking))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(perspective_taking)%>%
ggplot(aes(x = perspective_taking,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Perceived history of discrimination

Display code

summary(ds$history_discrimination)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       2       4       4       6       7      26 

Display code

ds %>% drop_na(history_discrimination)%>%
ggplot(aes(x = history_discrimination))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 10)+
  geom_textvline(label = "Mean = 4.00", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(history_discrimination, Country)%>%
ggplot(aes(x = history_discrimination))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 10)+
 labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(history_discrimination)%>%
ggplot(aes(x = history_discrimination,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Positive Affect

Display code

summary(ds$Event_positive_affect)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       1       2       4       6       7      32 

Display code

ds %>% drop_na(Event_positive_affect)%>%
ggplot(aes(x = Event_positive_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 4.00", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Positive affect score", 
       y = "Frequency", 
       title = "Event: Positive affect (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_positive_affect, Country)%>%
ggplot(aes(x = Event_positive_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Positive affect score", 
       y = "Frequency", 
       title = "Event: Positive affect by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_positive_affect)%>%
ggplot(aes(x = Event_positive_affect,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Positive affect score", 
       y = "Frequency", 
       title = "Event: Positive affect by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Negative Affect

Display code

summary(ds$Event_negative_affect)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     2.0     6.0     4.8     7.0     7.0      26 

Display code

ds %>% drop_na(Event_negative_affect)%>%
ggplot(aes(x = Event_negative_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Negative affect score", 
       y = "Frequency", 
       title = "Event: Negative affect (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_negative_affect, Country)%>%
ggplot(aes(x = Event_negative_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Negative affect score", 
       y = "Frequency", 
       title = "Event: Negative affect by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_negative_affect)%>%
ggplot(aes(x = Event_negative_affect,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Negative affect score", 
       y = "Frequency", 
       title = "Event: Negative affect by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Episodic recall

Display code

summary(ds$Event_episodic_recall)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       6       6       6       7       7      42 

Display code

ds %>% drop_na(Event_episodic_recall)%>%
ggplot(aes(x = Event_episodic_recall))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Episodic recall score", 
       y = "Frequency", 
       title = "Event: Episodic recall (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_episodic_recall, Country)%>%
ggplot(aes(x = Event_episodic_recall))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Episodic recall score", 
       y = "Frequency", 
       title = "Event: Episodic recall by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_episodic_recall)%>%
ggplot(aes(x = Event_episodic_recall,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Episodic recall score", 
       y = "Frequency", 
       title = "Event: Episodic recall by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Shared perception

Display code

summary(ds$Event_shared_perception)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      32 

Display code

ds %>% drop_na(Event_shared_perception)%>%
ggplot(aes(x = Event_shared_perception))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Shared perception score", 
       y = "Frequency", 
       title = "Event: Shared perception (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_shared_perception, Country)%>%
ggplot(aes(x = Event_shared_perception))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Shared perception score", 
       y = "Frequency", 
       title = "Event: Shared perception by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_shared_perception)%>%
ggplot(aes(x = Event_shared_perception,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Shared perception score", 
       y = "Frequency", 
       title = "Event: Shared perception by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Reflection

Display code

ds$Event_reflection <- ds$Event_event_reflection
summary(ds$Event_reflection)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       4       6       5       6       7      40 

Display code

ds %>% drop_na(Event_reflection)%>%
ggplot(aes(x = Event_reflection))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Reflection score", 
       y = "Frequency", 
       title = "Event: Reflection (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_reflection, Country)%>%
ggplot(aes(x = Event_reflection))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Reflection score", 
       y = "Frequency", 
       title = "Event: Reflection by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_reflection)%>%
ggplot(aes(x = Event_reflection,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Reflection score", 
       y = "Frequency", 
       title = "Event: Reflection by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Transformative for individual

Display code

summary(ds$Event_transformative_indiv)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       4       6       5       6       7      49 

Display code

ds %>% drop_na(Event_transformative_indiv)%>%
ggplot(aes(x = Event_transformative_indiv))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Transformative for individual score", 
       y = "Frequency", 
       title = "Event: Transformative for individual (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_transformative_indiv, Country)%>%
ggplot(aes(x = Event_transformative_indiv))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Transformative for individual score", 
       y = "Frequency", 
       title = "Event: Transformative for individual by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_transformative_indiv)%>%
ggplot(aes(x = Event_transformative_indiv,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Transformative for individual score", 
       y = "Frequency", 
       title = "Event: Transformative for individual by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Transformative for group

Display code

summary(ds$Event_transformative_group)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     4.0     5.0     4.9     6.0     7.0      28 

Display code

ds %>% drop_na(Event_transformative_group)%>%
ggplot(aes(x = Event_transformative_group))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Transformative for group score", 
       y = "Frequency", 
       title = "Event: Transformative for group (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_transformative_group, Country)%>%
ggplot(aes(x = Event_transformative_group))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Transformative for group score", 
       y = "Frequency", 
       title = "Event: Transformative for group by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_transformative_group)%>%
ggplot(aes(x = Event_transformative_group,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Transformative for group score", 
       y = "Frequency", 
       title = "Event: Transformative for group by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Imagistic event (sum of all subscales)

Display code

ds$Event_imagistic <- ds$Imagistic

summary(ds$Event_imagistic)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      7      32      36      35      40      49     102 

Display code

ds %>% drop_na(Event_imagistic)%>%
ggplot(aes(x = Event_imagistic))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 35.00", 
                 xintercept = 35.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Imagistic score", 
       y = "Frequency", 
       title = "Event: Imagistic (full sample)")+
  theme_bw()

Display code

ds %>% drop_na(Event_imagistic, Country)%>%
ggplot(aes(x = Event_imagistic))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Imagistic score", 
       y = "Frequency", 
       title = "Event: Imagistic by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code

ds %>% drop_na(Event_imagistic)%>%
ggplot(aes(x = Event_imagistic,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 35.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Imagistic score", 
       y = "Frequency", 
       title = "Event: Imagistic by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Section 3. Outcome: Social perception of Religious Freedom

Unconditional means model

Also called “varying intercept model with no predictors” (Gelman and Hill, 2016, Chapter 12). Allows intercepts to randomly vary across countries:

Display code

ds$SPRF <- ds$Sprf

## Varying intercept model with no predictors:
m00<- lmer(SPRF ~ 1 + (1 | Country), 
           data = ds)

summary(m00)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: SPRF ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 3068

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.726 -0.572  0.096  0.676  2.641 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0257   0.160   
 Residual             0.4502   0.671   
Number of obs: 1497, groups:  Country, 4

Fixed effects:
            Estimate Std. Error     df t value  Pr(>|t|)    
(Intercept)   5.4297     0.0822 3.0385      66 0.0000067 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random effects:

Variance for Intercept = 0.03. This is the variance of the means across level 1 categories (countries). Residual = 0.44. Variance explained by level 1 residuals (everything that’s not in level 1).

Display code

tab_model(m00)

	SPRF
Predictors	Estimates	CI	p
(Intercept)	5.43	5.27 – 5.59	<0.001
Random Effects
σ2	0.45
τ00 Country	0.03
ICC	0.05
N Country	4
Observations	1497
Marginal R2 / Conditional R2	0.000 / 0.054

We can see that ICC = 0.06. Lower ICC = low variance explained across groups. In this case, most of the variability is at individual-level (not group level). There is no significantly different patterns between countries.

Random intercept models

Also called “varying intercept model with individual-level predictors” (Gelman and Hill, 2016, Chapter 12).

Display code

## Create proper variable names:

ds$IG_Fusion <- (ds$IGF01+ds$IGF02+ds$IGF03)/3

ds$IG_Identification <- ds$igidentification

ds$OG_Bonds <- (ds$OGI01+ds$OGI02+ds$OGI03+
               ds$OGF01+ds$OGF02+ds$OGF03)/6

ds$Empathic_concern <- (ds$empathic_concern_01+ds$empathic_concern_02+
                        ds$empathic_concern_03)/3

ds$Perspective_taking <- (ds$perspective_taking_01+ds$perspective_taking_02+
                        ds$perspective_taking_03+ds$perspective_taking_03)/4

ds$Age <- ds$age

ds$Wealth_level <- as.factor(ifelse(ds$ses == "Lower middle", 1, 
                             ifelse(ds$ses == "Middle", 2,
                             ifelse(ds$ses == "Upper middle", 3,
                             ifelse(ds$ses == "Upper", 4, "")))))

ds$Female <- as.factor(ifelse(ds$gender == "Female", 1, 0))

ds$Married <- ds$married

ds$Perception_religious_freedom <- ds$SPRF

## Varying intercept models with individual-level predictors:
m01 <- lmer(Perception_religious_freedom~IG_Fusion+IG_Identification+OG_Bonds+Empathic_concern+
              Perspective_taking+Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m01)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Perception_religious_freedom ~ IG_Fusion + IG_Identification +  
    OG_Bonds + Empathic_concern + Perspective_taking + Age +  
    Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 2526

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-5.087 -0.593  0.033  0.677  4.070 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0192   0.139   
 Residual             0.3797   0.616   
Number of obs: 1315, groups:  Country, 4

Fixed effects:
                      Estimate  Std. Error          df t value
(Intercept)           2.791389    0.186736  126.440741   14.95
IG_Fusion             0.076776    0.022317 1301.703278    3.44
IG_Identification     0.124264    0.024302 1299.630921    5.11
OG_Bonds              0.017347    0.011905 1301.910330    1.46
Empathic_concern      0.038961    0.018456 1301.259994    2.11
Perspective_taking    0.174826    0.020026 1301.014809    8.73
Age                  -0.000787    0.001384 1301.732453   -0.57
Female1               0.039063    0.035127 1301.240598    1.11
MarriedOther         -0.015903    0.081643 1300.212099   -0.19
MarriedUnmarried      0.088790    0.040799 1300.996291    2.18
Wealth_level2         0.164497    0.044162 1201.020084    3.72
Wealth_level3         0.138814    0.063917 1244.155630    2.17
Wealth_level4        -0.035545    0.130529 1299.892380   -0.27
                               Pr(>|t|)    
(Intercept)        < 0.0000000000000002 ***
IG_Fusion                        0.0006 ***
IG_Identification            0.00000036 ***
OG_Bonds                         0.1453    
Empathic_concern                 0.0350 *  
Perspective_taking < 0.0000000000000002 ***
Age                              0.5694    
Female1                          0.2663    
MarriedOther                     0.8456    
MarriedUnmarried                 0.0297 *  
Wealth_level2                    0.0002 ***
Wealth_level3                    0.0301 *  
Wealth_level4                    0.7854    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m01)

	Perception_religious_freedom
Predictors	Estimates	CI	p
(Intercept)	2.79	2.43 – 3.16	<0.001
IG Fusion	0.08	0.03 – 0.12	0.001
IG Identification	0.12	0.08 – 0.17	<0.001
OG Bonds	0.02	-0.01 – 0.04	0.145
Empathic concern	0.04	0.00 – 0.08	0.035
Perspective taking	0.17	0.14 – 0.21	<0.001
Age	-0.00	-0.00 – 0.00	0.569
Female [1]	0.04	-0.03 – 0.11	0.266
Married [Other]	-0.02	-0.18 – 0.14	0.846
Married [Unmarried]	0.09	0.01 – 0.17	0.030
Wealth level [2]	0.16	0.08 – 0.25	<0.001
Wealth level [3]	0.14	0.01 – 0.26	0.030
Wealth level [4]	-0.04	-0.29 – 0.22	0.785
Random Effects
σ2	0.38
τ00 Country	0.02
ICC	0.05
N Country	4
Observations	1315
Marginal R2 / Conditional R2	0.170 / 0.210

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code

## Change class of all models so we can use stargazer():
class(m00) <- "lmerMod"
class(m01) <- "lmerMod"

## Tabulated results:
stargazer(m00, m01,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	SPRF	Perception_religious_freedom
	(1)	(2)
IG_Fusion		0.077***
		(0.022)
IG_Identification		0.120***
		(0.024)
OG_Bonds		0.017
		(0.012)
Empathic_concern		0.039*
		(0.018)
Perspective_taking		0.170***
		(0.020)
Age		-0.001
		(0.001)
Female1		0.039
		(0.035)
MarriedOther		-0.016
		(0.082)
MarriedUnmarried		0.089*
		(0.041)
Wealth_level2		0.160***
		(0.044)
Wealth_level3		0.140*
		(0.064)
Wealth_level4		-0.036
		(0.130)
Constant	5.400***	2.800***
	(0.082)	(0.190)
Observations	1,497	1,315
Log Likelihood	-1,534.000	-1,263.000
Akaike Inf. Crit.	3,074.000	2,556.000
Bayesian Inf. Crit.	3,090.000	2,634.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Random intercept models: Imagistic predictors

Display code

## Varying intercept models with individual-level predictors:
m02 <- lmer(Perception_religious_freedom~Event_shared_perception+Event_episodic_recall+
              Event_reflection+Event_positive_affect+Event_negative_affect+
              Event_transformative_indiv+Event_transformative_group+
              Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m02)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
Perception_religious_freedom ~ Event_shared_perception + Event_episodic_recall +  
    Event_reflection + Event_positive_affect + Event_negative_affect +  
    Event_transformative_indiv + Event_transformative_group +  
    Age + Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 2640

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.546 -0.647  0.025  0.683  3.091 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0359   0.189   
 Residual             0.3782   0.615   
Number of obs: 1369, groups:  Country, 4

Fixed effects:
                               Estimate   Std. Error           df t value
(Intercept)                   3.6878965    0.1683883   26.5844943   21.90
Event_shared_perception       0.0182817    0.0149981 1353.5600531    1.22
Event_episodic_recall         0.1463361    0.0180858 1352.8118832    8.09
Event_reflection              0.0071473    0.0147259 1343.1801914    0.49
Event_positive_affect         0.0127756    0.0098358 1353.0070788    1.30
Event_negative_affect         0.0343717    0.0094891 1353.3489940    3.62
Event_transformative_indiv    0.0448096    0.0152013 1353.6154566    2.95
Event_transformative_group    0.0369059    0.0139939 1352.2634171    2.64
Age                           0.0000562    0.0013581 1352.6211111    0.04
Female1                       0.0262406    0.0342730 1353.8540285    0.77
MarriedOther                 -0.0436935    0.0815275 1353.6452664   -0.54
MarriedUnmarried              0.0826492    0.0399443 1351.8882536    2.07
Wealth_level2                 0.1564298    0.0436726 1318.8752569    3.58
Wealth_level3                 0.0181246    0.0628422 1336.7050366    0.29
Wealth_level4                -0.1740043    0.1303406 1353.9648284   -1.33
                                       Pr(>|t|)    
(Intercept)                < 0.0000000000000002 ***
Event_shared_perception                 0.22308    
Event_episodic_recall        0.0000000000000013 ***
Event_reflection                        0.62751    
Event_positive_affect                   0.19421    
Event_negative_affect                   0.00030 ***
Event_transformative_indiv              0.00326 ** 
Event_transformative_group              0.00845 ** 
Age                                     0.96703    
Female1                                 0.44403    
MarriedOther                            0.59209    
MarriedUnmarried                        0.03873 *  
Wealth_level2                           0.00035 ***
Wealth_level3                           0.77307    
Wealth_level4                           0.18210    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m02)

	Perception_religious_freedom
Predictors	Estimates	CI	p
(Intercept)	3.69	3.36 – 4.02	<0.001
Event shared perception	0.02	-0.01 – 0.05	0.223
Event episodic recall	0.15	0.11 – 0.18	<0.001
Event reflection	0.01	-0.02 – 0.04	0.628
Event positive affect	0.01	-0.01 – 0.03	0.194
Event negative affect	0.03	0.02 – 0.05	<0.001
Event transformative indiv	0.04	0.01 – 0.07	0.003
Event transformative group	0.04	0.01 – 0.06	0.008
Age	0.00	-0.00 – 0.00	0.967
Female [1]	0.03	-0.04 – 0.09	0.444
Married [Other]	-0.04	-0.20 – 0.12	0.592
Married [Unmarried]	0.08	0.00 – 0.16	0.039
Wealth level [2]	0.16	0.07 – 0.24	<0.001
Wealth level [3]	0.02	-0.11 – 0.14	0.773
Wealth level [4]	-0.17	-0.43 – 0.08	0.182
Random Effects
σ2	0.38
τ00 Country	0.04
ICC	0.09
N Country	4
Observations	1369
Marginal R2 / Conditional R2	0.145 / 0.219

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code

## Change class of all models so we can use stargazer():
class(m00) <- "lmerMod"
class(m01) <- "lmerMod"
class(m02) <- "lmerMod"

## Tabulated results:
stargazer(m00, m01, m02,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	SPRF	Perception_religious_freedom
	(1)	(2)	(3)
IG_Fusion		0.077***
		(0.022)
IG_Identification		0.120***
		(0.024)
OG_Bonds		0.017
		(0.012)
Empathic_concern		0.039*
		(0.018)
Perspective_taking		0.170***
		(0.020)
Event_shared_perception			0.018
			(0.015)
Event_episodic_recall			0.150***
			(0.018)
Event_reflection			0.007
			(0.015)
Event_positive_affect			0.013
			(0.010)
Event_negative_affect			0.034***
			(0.009)
Event_transformative_indiv			0.045**
			(0.015)
Event_transformative_group			0.037**
			(0.014)
Age		-0.001	0.0001
		(0.001)	(0.001)
Female1		0.039	0.026
		(0.035)	(0.034)
MarriedOther		-0.016	-0.044
		(0.082)	(0.082)
MarriedUnmarried		0.089*	0.083*
		(0.041)	(0.040)
Wealth_level2		0.160***	0.160***
		(0.044)	(0.044)
Wealth_level3		0.140*	0.018
		(0.064)	(0.063)
Wealth_level4		-0.036	-0.170
		(0.130)	(0.130)
Constant	5.400***	2.800***	3.700***
	(0.082)	(0.190)	(0.170)
Observations	1,497	1,315	1,369
Log Likelihood	-1,534.000	-1,263.000	-1,320.000
Akaike Inf. Crit.	3,074.000	2,556.000	2,673.000
Bayesian Inf. Crit.	3,090.000	2,634.000	2,762.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Perception of Religious Freedom

Display code

summary(ds$SPRF)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      2       5       6       5       6       7      45 

Display code

ggplot(data = ds, 
       aes(x = SPRF)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.0", 
                 xintercept = 5.0, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "SPRF score", 
       title = "Social Perception of Religiour Freedom")+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(x = SPRF)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "SPRF score", 
       title = "Social Perception of Religiour Freedom")+
  facet_wrap( ~ Country, nrow = 2) +
  theme_bw()

Display code

df01 <- ds %>% drop_na(SPRF)

tbl01 <- aggregate(df01$SPRF, 
                    by=list(df01$Country),
                    FUN=mean)
tbl01$Country <- tbl01$Group.1
tbl01$SPRF <- tbl01$x
tbl01 <- tbl01[, 3:4]
tbl01

   Country SPRF
1   Gambia  5.6
2 Pakistan  5.2
3 Tanzania  5.4
4   Uganda  5.6

Display code

ggplot(data = df01, 
       aes(x = SPRF, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 4.70", 
                 xintercept = 4.70, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "SPRF score", 
       title = "Social Perceptions of Relgious Freedom")+
  theme_bw()

Section 3. Outcome: Experience of Religious Freedom

Unconditional means model

Also called “varying intercept model with no predictors” (Gelman and Hill, 2016, Chapter 12). Allows intercepts to randomly vary across countries:

Display code

ds$Exp_religious_freedom <- ds$Exp_religious_freedom

## Varying intercept model with no predictors:
m10<- lmer(Exp_religious_freedom ~ 1 + (1 | Country), 
           data = ds)

summary(m10)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Exp_religious_freedom ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 5077

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.924 -0.515  0.369  0.705  1.210 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.129    0.36    
 Residual             1.743    1.32    
Number of obs: 1492, groups:  Country, 4

Fixed effects:
            Estimate Std. Error    df t value Pr(>|t|)    
(Intercept)    5.841      0.183 3.037    31.9 0.000062 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random effects:

Variance for Intercept = 0.141. This is the variance of the means across level 1 categories (countries). Residual = 1.7. Variance explained by level 1 residuals (everything that’s not in level 1).

Display code

tab_model(m10)

	Exp_religious_freedom
Predictors	Estimates	CI	p
(Intercept)	5.84	5.48 – 6.20	<0.001
Random Effects
σ2	1.74
τ00 Country	0.13
ICC	0.07
N Country	4
Observations	1492
Marginal R2 / Conditional R2	0.000 / 0.069

We can see that ICC = 0.08. Lower ICC = low variance explained across groups. In this case, most of the variability is at individual-level (not group level). There is no significantly different patterns between countries.

Random intercept models

Also called “varying intercept model with individual-level predictors” (Gelman and Hill, 2016, Chapter 12).

Display code

## Varying intercept models with individual-level predictors:
m11 <- lmer(Exp_religious_freedom~IG_Fusion+IG_Identification+OG_Bonds+Empathic_concern+
              Perspective_taking+Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m11)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Exp_religious_freedom ~ IG_Fusion + IG_Identification + OG_Bonds +  
    Empathic_concern + Perspective_taking + Age + Female + Married +  
    Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 4251

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.592 -0.524  0.233  0.696  2.405 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0852   0.292   
 Residual             1.4240   1.193   
Number of obs: 1316, groups:  Country, 4

Fixed effects:
                     Estimate Std. Error         df t value
(Intercept)           4.50051    0.36592   98.91990   12.30
IG_Fusion             0.04887    0.04321 1302.16380    1.13
IG_Identification     0.11225    0.04709 1300.48983    2.38
OG_Bonds             -0.20987    0.02305 1302.97689   -9.10
Empathic_concern      0.17408    0.03574 1302.85029    4.87
Perspective_taking   -0.00743    0.03884 1302.50507   -0.19
Age                   0.00260    0.00267 1302.57372    0.97
Female1               0.22350    0.06807 1302.67705    3.28
MarriedOther          0.22189    0.15809 1302.28377    1.40
MarriedUnmarried      0.04714    0.07906 1301.68725    0.60
Wealth_level2        -0.03748    0.08541 1235.25598   -0.44
Wealth_level3        -0.34497    0.12377 1264.82102   -2.79
Wealth_level4        -0.49288    0.25279 1302.08011   -1.95
                               Pr(>|t|)    
(Intercept)        < 0.0000000000000002 ***
IG_Fusion                        0.2583    
IG_Identification                0.0173 *  
OG_Bonds           < 0.0000000000000002 ***
Empathic_concern              0.0000013 ***
Perspective_taking               0.8484    
Age                              0.3303    
Female1                          0.0011 ** 
MarriedOther                     0.1607    
MarriedUnmarried                 0.5512    
Wealth_level2                    0.6609    
Wealth_level3                    0.0054 ** 
Wealth_level4                    0.0514 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m11)

	Exp_religious_freedom
Predictors	Estimates	CI	p
(Intercept)	4.50	3.78 – 5.22	<0.001
IG Fusion	0.05	-0.04 – 0.13	0.258
IG Identification	0.11	0.02 – 0.20	0.017
OG Bonds	-0.21	-0.26 – -0.16	<0.001
Empathic concern	0.17	0.10 – 0.24	<0.001
Perspective taking	-0.01	-0.08 – 0.07	0.848
Age	0.00	-0.00 – 0.01	0.330
Female [1]	0.22	0.09 – 0.36	0.001
Married [Other]	0.22	-0.09 – 0.53	0.161
Married [Unmarried]	0.05	-0.11 – 0.20	0.551
Wealth level [2]	-0.04	-0.21 – 0.13	0.661
Wealth level [3]	-0.34	-0.59 – -0.10	0.005
Wealth level [4]	-0.49	-0.99 – 0.00	0.051
Random Effects
σ2	1.42
τ00 Country	0.09
ICC	0.06
N Country	4
Observations	1316
Marginal R2 / Conditional R2	0.122 / 0.172

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code

## Change class of all models so we can use stargazer():
class(m10) <- "lmerMod"
class(m11) <- "lmerMod"

## Tabulated results:
stargazer(m10, m11,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Exp_religious_freedom
	(1)	(2)
IG_Fusion		0.049
		(0.043)
IG_Identification		0.110*
		(0.047)
OG_Bonds		-0.210***
		(0.023)
Empathic_concern		0.170***
		(0.036)
Perspective_taking		-0.007
		(0.039)
Age		0.003
		(0.003)
Female1		0.220**
		(0.068)
MarriedOther		0.220
		(0.160)
MarriedUnmarried		0.047
		(0.079)
Wealth_level2		-0.037
		(0.085)
Wealth_level3		-0.340**
		(0.120)
Wealth_level4		-0.490
		(0.250)
Constant	5.800***	4.500***
	(0.180)	(0.370)
Observations	1,492	1,316
Log Likelihood	-2,538.000	-2,125.000
Akaike Inf. Crit.	5,083.000	4,281.000
Bayesian Inf. Crit.	5,099.000	4,359.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Random intercept models: Imagistic predictors

Display code

## Varying intercept models with individual-level predictors:
m12 <- lmer(Exp_religious_freedom~Event_shared_perception+Event_episodic_recall+
              Event_reflection+Event_positive_affect+Event_negative_affect+
              Event_transformative_indiv+Event_transformative_group+
              Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m12)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
Exp_religious_freedom ~ Event_shared_perception + Event_episodic_recall +  
    Event_reflection + Event_positive_affect + Event_negative_affect +  
    Event_transformative_indiv + Event_transformative_group +  
    Age + Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 4554

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.431 -0.583  0.280  0.722  1.891 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.126    0.355   
 Residual             1.571    1.253   
Number of obs: 1365, groups:  Country, 4

Fixed effects:
                             Estimate Std. Error         df t value
(Intercept)                   5.99486    0.33420   33.29321   17.94
Event_shared_perception       0.05464    0.03048 1349.92882    1.79
Event_episodic_recall         0.23886    0.03705 1349.09559    6.45
Event_reflection             -0.12720    0.03041 1329.62089   -4.18
Event_positive_affect        -0.09007    0.02018 1347.71554   -4.46
Event_negative_affect        -0.14057    0.01943 1349.59958   -7.23
Event_transformative_indiv    0.00229    0.03112 1349.11898    0.07
Event_transformative_group   -0.08194    0.02868 1348.51282   -2.86
Age                           0.00420    0.00276 1348.89116    1.52
Female1                       0.22681    0.06997 1349.95972    3.24
MarriedOther                  0.32713    0.16619 1349.86561    1.97
MarriedUnmarried              0.15721    0.08155 1348.03301    1.93
Wealth_level2                -0.14879    0.08869 1302.56527   -1.68
Wealth_level3                -0.45868    0.12821 1325.58354   -3.58
Wealth_level4                -0.68178    0.26560 1349.76810   -2.57
                                       Pr(>|t|)    
(Intercept)                < 0.0000000000000002 ***
Event_shared_perception                 0.07325 .  
Event_episodic_recall          0.00000000015806 ***
Event_reflection               0.00003074047694 ***
Event_positive_affect          0.00000871652381 ***
Event_negative_affect          0.00000000000078 ***
Event_transformative_indiv              0.94143    
Event_transformative_group              0.00434 ** 
Age                                     0.12849    
Female1                                 0.00122 ** 
MarriedOther                            0.04923 *  
MarriedUnmarried                        0.05409 .  
Wealth_level2                           0.09367 .  
Wealth_level3                           0.00036 ***
Wealth_level4                           0.01037 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m12)

	Exp_religious_freedom
Predictors	Estimates	CI	p
(Intercept)	5.99	5.34 – 6.65	<0.001
Event shared perception	0.05	-0.01 – 0.11	0.073
Event episodic recall	0.24	0.17 – 0.31	<0.001
Event reflection	-0.13	-0.19 – -0.07	<0.001
Event positive affect	-0.09	-0.13 – -0.05	<0.001
Event negative affect	-0.14	-0.18 – -0.10	<0.001
Event transformative indiv	0.00	-0.06 – 0.06	0.941
Event transformative group	-0.08	-0.14 – -0.03	0.004
Age	0.00	-0.00 – 0.01	0.128
Female [1]	0.23	0.09 – 0.36	0.001
Married [Other]	0.33	0.00 – 0.65	0.049
Married [Unmarried]	0.16	-0.00 – 0.32	0.054
Wealth level [2]	-0.15	-0.32 – 0.03	0.094
Wealth level [3]	-0.46	-0.71 – -0.21	<0.001
Wealth level [4]	-0.68	-1.20 – -0.16	0.010
Random Effects
σ2	1.57
τ00 Country	0.13
ICC	0.07
N Country	4
Observations	1365
Marginal R2 / Conditional R2	0.095 / 0.162

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code

## Change class of all models so we can use stargazer():
class(m10) <- "lmerMod"
class(m11) <- "lmerMod"
class(m12) <- "lmerMod"

## Tabulated results:
stargazer(m10, m11, m12,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Exp_religious_freedom
	(1)	(2)	(3)
IG_Fusion		0.049
		(0.043)
IG_Identification		0.110*
		(0.047)
OG_Bonds		-0.210***
		(0.023)
Empathic_concern		0.170***
		(0.036)
Perspective_taking		-0.007
		(0.039)
Event_shared_perception			0.055
			(0.030)
Event_episodic_recall			0.240***
			(0.037)
Event_reflection			-0.130***
			(0.030)
Event_positive_affect			-0.090***
			(0.020)
Event_negative_affect			-0.140***
			(0.019)
Event_transformative_indiv			0.002
			(0.031)
Event_transformative_group			-0.082**
			(0.029)
Age		0.003	0.004
		(0.003)	(0.003)
Female1		0.220**	0.230**
		(0.068)	(0.070)
MarriedOther		0.220	0.330*
		(0.160)	(0.170)
MarriedUnmarried		0.047	0.160
		(0.079)	(0.082)
Wealth_level2		-0.037	-0.150
		(0.085)	(0.089)
Wealth_level3		-0.340**	-0.460***
		(0.120)	(0.130)
Wealth_level4		-0.490	-0.680*
		(0.250)	(0.270)
Constant	5.800***	4.500***	6.000***
	(0.180)	(0.370)	(0.330)
Observations	1,492	1,316	1,365
Log Likelihood	-2,538.000	-2,125.000	-2,277.000
Akaike Inf. Crit.	5,083.000	4,281.000	4,588.000
Bayesian Inf. Crit.	5,099.000	4,359.000	4,677.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Experience of Religious Freedom

Display code

df01 <- ds%>% drop_na(Exp_religious_freedom)
summary(df01$Exp_religious_freedom)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     5.1     6.3     5.8     7.0     7.0 

Display code

ggplot(data = df01, 
       aes(x = Exp_religious_freedom)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.8", 
                 xintercept = 5.8, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "Exp_religious_freedom score", 
       title = "Experience of Relgious Freedom")+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(x = Exp_religious_freedom)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "Exp_religious_freedom score", 
       title = "Experience of Religious Freedom")+
  facet_wrap( ~ Country, nrow = 2) +
  theme_bw()

Display code

df01 <- ds %>% drop_na(Exp_religious_freedom)

tbl01 <- aggregate(df01$Exp_religious_freedom, 
                    by=list(df01$Country),
                    FUN=mean)
tbl01$Country <- tbl01$Group.1
tbl01$Exp_religious_freedom <- tbl01$x
tbl01 <- tbl01[, 3:4]
tbl01

   Country Exp_religious_freedom
1   Gambia                   6.1
2 Pakistan                   5.4
3 Tanzania                   5.7
4   Uganda                   6.2

Display code

ggplot(data = df01, 
       aes(x = Exp_religious_freedom, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.8", 
                 xintercept = 5.8, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "Exp_religious_freedom score", 
       title = "Experience of Relgious Freedom")+
  theme_bw()

Section X: Coefficient Plots

Outcome: Social Perception of Religious Freedom

Display code

mp01 <- modelplot(m01,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Social perception of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp01

Display code

mp02 <- modelplot(m02,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Social perception of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp02

Display code

mp01a <- mp01 + theme(legend.position = "none")
ol2 <- mp01a / mp02
ol2

Outcome: Experience of Religious Freedom

Display code

mp11 <- modelplot(m11,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Experience of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp11

Display code

mp12 <- modelplot(m12,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Experience of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp12

Display code

mp11a <- mp11 + theme(legend.position = "none")
ol2a <- mp11a / mp12
ol2a